REM ------ Q2Ja --------- Julia Set of a rational function ----------------
REM
REM              Q2(z)=z/(z^2-az+a),   z=x+iy, a=ar+i*ai
REM  
REM     与えられた四角形の領域を600×600の格子点に分割し，各点を初期値
REM    として反復，１つの吸引的周期軌道に収束するとき，初期値に同じ色を
REM    つける。吸引的周期軌道が２つあることが前提。
REM    塗り絵は2色：
REM          深緑   ．．．1つ目の吸引的軌道に収束  jo=1
REM          赤(or 黄) ．．．2つ目の        〃     jo=2
REM          オリーブ色(13)  ．．．収束判定不明    jo=0
REM        各色の境界がジュリア集合　(注）カラーは時々変える
REM
REM   2つの吸引的周期軌道の周期および周期点のうち１点を入力する（計2点）。
REM                q=qr(k)+i*qi(k),  k=1,2
REM           収束判定は  |Zn - ｑ|<0.001
REM     
REM           it(1),it(2): 各色のドット数  
REM           it(0): 収束判定不明のドット数
REM-----------------------------------  7.2 '19  joe  ---------------------   
DIM co(0 TO 2),it(0 TO 2),qr(2),qi(2),per(2)
READ co(1),co(2)
DATA 10,4
REM DATA 6,10
LINE INPUT PROMPT"テスト(t)  or  本番(h)？ t/h=":tq$
IF tq$="H" OR tq$="h" THEN GOTO 10
READ txl,txu,tyl,tyu,n1
DATA .3,.5,.6,.8,200
GOTO 20
10 LET txl=.1
   LET txu=.9
   LET tyl=.1
   LET tyu=.9
   LET n1=600
20 REM ---------
   FOR i=1 TO 2
      LET it(i)=0
   NEXT i
   LET EPS=.000001
   INPUT PROMPT "パラメタaのinput (ar,ai)=":ar,ai
   print"複素平面上の四角形のinput"
   INPUT PROMPT "xの範囲 xl,xu=":xl,xu
   INPUT PROMPT "yの範囲 yl,yu=":yl,yu
   REM ---- 画像パラメタ設定  ----
   LET wx=xu-xl
   LET wy=yu-yl
   LET twx=txu-txl
   LET twy=tyu-tyl
   LET a1=twx/wx
   LET a2=twy/wy
   LET b1=(xu*txl-xl*txu)/wx
   LET b2=(yu*tyl-yl*tyu)/wy
   SET AREA COLOR 13
   PLOT AREA:txl,tyl;txu,tyl;txu,tyu;txl,tyu;txl,tyl
   SET LINE COLOR 7
   IF xl*xu<0 THEN PLOT LINES:b1,tyl-.02;b1,tyu+.02
   IF yu*yl<0 THEN PLOT LINES:txl-.02,b2;txu+.02,b2
   PRINT"周期軌道の周期と周期点をinputせよ"
   INPUT PROMPT"最初のものの周期 per(1)=":per(1)
   INPUT PROMPT"周期点(qr(1),qi(1))=":qr(1),qi(1)
   INPUT PROMPT"2番目のものの周期 per(2)=":per(2)
   INPUT PROMPT"周期点(qr(2),qi(2))=":qr(2),qi(2)
   INPUT PROMPT"iteration の回数(tm 100以上)=":tm
   REM ------------------------------------------------------
   PLOT TEXT, AT .04,.96:"Q2-Julia Set,  a="
   PLOT TEXT, AT .27,.96,USING"##.##### ##.#####":ar,ai
   PLOT TEXT, AT  .6,.96:"反復回数 tm="
   PLOT TEXT, AT .76,.96,USING"####":tm
   PLOT TEXT, AT.04,.92:"複素平面:"
   PLOT TEXT, AT .22,.92,USING"#####.##### #####.#####":xl,xu
   PLOT TEXT, AT .56,.92,USING"#####.##### #####.#####":yl,yu
   PLOT TEXT, AT .04,.06:"周期, 周期点="
   PLOT TEXT, AT .2,.06,USING"##  (##.####### ##.#######)":per(1),qr(1),qi(1)
   PLOT TEXT, AT .55,.06,USING"##  (##.####### ##.#######)":per(2),qr(2),qi(2) 
   REM
   REM
   REM     計算開始
   REM
   LET IS=0
   LET dx=wx/n1
   LET dy=wy/n1
   SET POINT STYLE 1
   FOR i=0 TO n1
      FOR j=0 TO n1
         LET x=xl+i*dx
         LET tx=a1*x+b1
         LET y=yl+j*dy
         LET ty=a2*y+b2
         REM
         LET jo=0
         FOR t=1 TO tm 
            LET d1=x*x-y*y-ar*x+ai*y+ar
            LET d2=2*x*y-ai*x-ar*y+ai
            LET d=d1*d1+d2*d2
            IF d<EPS THEN LET IS=IS+1
            LET x1=(x*d1+y*d2)/d
            LET y1=(y*d1-x*d2)/d
            IF t<60 THEN GOTO 40 ELSE GOTO 50
40          LET x=x1
            LET y=y1
            GOTO 70
50          FOR k=1 TO 2
               LET e=SQR((x1-qr(k))^2+(y1-qi(k))^2)
               IF e<0.001 THEN GOTO 55 ELSE GOTO 60
55             LET jo=k
               GOTO 80
60          NEXT k
            LET x=x1
            LET y=y1
70       NEXT t
80       REM  ---- 塗り絵  -----
         LET it(jo)=it(jo)+1
         IF jo=0 THEN GOTO 85
         SET POINT COLOR co(jo)   
         PLOT  POINTS: tx,ty
85    NEXT j
   NEXT i
   REM  ----------------------------------------------------
   PRINT"分母ほとんどゼロの回数=";IS
   REM
   REM ---------------------------
   INPUT PROMPT"ファイル名，コメント書くか? y/n ":q$
   IF q$="n" OR q$="n" THEN GOTO 90
   PRINT "ファイル名などinputせよ(全角20文字以下) (例）Q2ja-01,Q2ja-02,.."
   LINE INPUT PROMPT "ファイル名など:":na$
   PLOT TEXT, AT .02,.02:na$
90 PLOT TEXT,AT .45,.02:"co(1),co(2),ｵﾘｰﾌﾞ:"
   REM 90 PLOT TEXT,AT .42,.02:"銀色,えび茶,ｵﾘｰﾌﾞ:"
   PLOT TEXT, AT .68,.02,USING"###### ###### ######":it(1),it(2),it(0)
END